capture log close analysis
log using mml-analysis, name(analysis) replace text

//  program:    mml-analysis.do
//  task:		replication of Hall medical marijuana study
// 	input:		ndush.dta
//	output:		mml-results.txt
//  project:    medical marijuana laws
//  author:     sam harper \ 02dec2011

//  #0
//  program setup

version 11
set linesize 80
clear all
macro drop _all

set matsize 1100
set more off

// #1
// bring in prevalence data from National Survey of Drug Use & Health

use mml-data.dta, clear

// #2
// replication of Wall et al. estimates

* limit to same years as Wall et al.
keep if year>=2002 & year<=2007

* overall prevalence
sum pmu if age==2

* stratified by ever treated (extended definition)
tab treatex if age==2, sum(pmu)

* Wall et al. Table 1
table year0 mmlf if age==2, c(mean pmu)
table year0 mmlf if age==2, c(mean pru)

* test for interaction of treatment status with year fixed effects
reg pmu year0##mmlf if age==2
testparm i(1/5).year0#1.mmlf

reg pru year0##mmlf if age==2
testparm i(1/5).year0#1.mmlf

* Wall et al. random effects regression estimates
xtreg pmu year0 i.wall if age==2, i(state) re
margins wall

xtreg pru year0 i.wall if age==2, i(state) re
margins wall

* Fixed effects for year vs. diff-in-diffs
foreach var of varlist pmu pru {
	
	* no state fixed effects
	reg `var' i.year0 i.mmlf if age==2, cformat(%4.3f) 
	scalar wallrep=_b[1.mmlf] // save treatment effect from Wall model
	
	* difference-in-differences
	xtreg `var' i.year0 i.mmlf if age==2, i(state) fe cluster(state) cformat(%4.3f)
	test 1.mmlf=wallrep
	
	* difference-in-differences with restricted control group
	xtreg `var' i.year0 i.mmlf if age==2 & treatex==1, i(state) fe cluster(state) cformat(%4.3f)
	}

* t-tests for comparisons, written to file
statsby nstates=r(N_2) wmml=r(mu_2) womml=r(mu_1) tstat=r(t) ///
 pval=r(p), by(year) saving(wallT1a, replace): ttest pmu if age==2, by(mmlf)
use wallT1a, clear
outsheet using wallT1a.csv, replace

use mml-data.dta, clear
statsby nstates=r(N_2) wmml=r(mu_2) womml=r(mu_1) tstat=r(t) ///
 pval=r(p), by(year) saving(wallT1b, replace): ttest pru if age==2, by(mmlf)
use wallT1b, clear
outsheet using wallT1b.csv, replace

// #2
// Replication of regression and difference-in-differences estimates
use mml-data.dta, clear

* set up matrix for collecting test results
	matrix pmu_walldd_comp = J(4,4,-99)
	matrix pru_walldd_comp = J(4,4,-99)
	matrix colnames pmu_walldd_comp= age wall dd pval_wall_dd
	matrix colnames pru_walldd_comp= age wall dd pval_wall_dd


* for different age groups and two different outcomes
foreach var of varlist pmu pru {
qui levelsof age, local(levels)
foreach a of local levels {
    di " "
    di "-> Age group = `: label (age) `a''"
    
    * no state fixed effects
	reg `var' i.year0 i.mmlf if age==`a', cformat(%4.3f) 
	scalar ewall=_b[1.mmlf] // save treatment effect from Wall model
	outreg2 using mml-`var', sideway stats(coef ci) paren(ci) dec(2) cdec(1) noaster ///
		addtext(Year fixed effects, Yes, State fixed effects, No, Age group, `a')
	
	matrix `var'_walldd_comp[`a',2]=ewall
		
	* difference-in-differeces model
	xtreg `var' i.year0 i.mmlf if age==`a', i(state) fe cluster(state) cformat(%4.3f)
	scalar firsty=_b[1.mmlf] // save fixed effect estimate with FE for each year
	outreg2 using mml-`var', sideway stats(coef ci) paren(ci) dec(2) cdec(1) noaster ///
		addtext(Year fixed effects, Yes, State fixed effects, Yes, Age group, `a')
	
	matrix `var'_walldd_comp[`a',3]=firsty
	
	* test of whether treatment effect differs between fixed effect model
	* and random effect model as in Wall et el.
	di "Test of DD vs. estimate of Wall et al.: "
	test 1.mmlf=ewall
	scalar `var'femmlf`a'=_b[1.mmlf]
	
	*place test results in matrix
	matrix rownames `var'_walldd_comp = `var'
	matrix `var'_walldd_comp[`a',4]=r(p)
	matrix `var'_walldd_comp[`a',1]=`a'
		}
	}

di "Summary of tests of fixed vs. random effect models"
di "Prevalence of monthly use: "
mat l pmu_walldd_comp, format(%6.4g)
di " "
di "Perceived riskiness of monthly use: "
mat l pru_walldd_comp, format(%6.4g)

// #3
// incorporate error in outcome measurement
* define program to draw randomly from the measured outcome
* and run the diff-in-diffs regression
program drop _all
program mmlsim, rclassversion 11
args yvar agen `yvar'sim = rnormal(`yvar',`yvar'se)
xtreg `yvar'sim i.mmlf i.year0 if age==`a', i(state) fe cluster(state) cformat(%4.3f)
drop `yvar'sim
end

* simulate and grab coefficients
use mml-data, clear
foreach var of varlist pmu pru {
	qui levelsof age, local(levels)
	foreach i of local levels {
    di " "
    di "-> Age group = `: label (age) `i''"
		simulate coef`var'`i'=_b[1.mmlf], nodots seed(72903) ///
		reps(5000) saving(mmlsim`var'`i', replace): mmlsim `var' `i'
		sum coef`var'`i'
		centile coef`var'`i', centile(2.5, 97.5)
		use mml-data, clear
		}
	}


// #4
// sensitivity analysis for assuming law affects use in second
// rather than first year of observed data


* set up matrix for collecting test results
	matrix pmu_year_comp = J(4,4,-99)
	matrix pru_year_comp = J(4,4,-99)
	matrix colnames pmu_year_comp= age mmlf mmll pval
	matrix colnames pru_year_comp= age mmlf mmll pval

* for different age groups and two different outcomes
foreach var of varlist pmu pru {
qui levelsof age, local(levels)
foreach a of local levels {
    di " "
    di "-> Age group = `: label (age) `a''"

	qui xtreg `var' i.year0 i.mmlf if age==`a', i(state) fe cluster(state) cformat(%4.3f)
	scalar `var'femmlf`a'=_b[1.mmlf]
	xtreg `var' i.year0 i.mmll if age==`a', i(state) fe cluster(state) cformat(%4.3f)
	outreg2 using mml-`var', sideway stats(coef ci) paren(ci) dec(2) cdec(1) noaster ///
		addtext(Year fixed effects, Yes, State fixed effects, Yes, Age group, `a')	
	
	* test of whether treatment effect differs between fixed effect model
	* with first (mmlf) vs. last (mmll) year of implementation
	di "Test of first vs. last year of implementation: "
	test 1.mmll=`var'femmlf`a'
	
	* matrix for collecting test results
	matrix rownames `var'_year_comp = `var'
	matrix `var'_year_comp[`a',1]=`a'
	matrix `var'_year_comp[`a',2]=`var'femmlf`a'
	scalar lasty=_b[1.mmll]
	matrix `var'_year_comp[`a',3]=lasty
	matrix `var'_year_comp[`a',4]=r(p)
		}
	}
di "Summary of tests of first vs. last year of implementation"
di "Prevalence of monthly use: "
mat l pmu_year_comp, format(%6.4g)
di " "
di "Perceived riskiness of monthly use: "
mat l pru_year_comp, format(%6.4g)

// #5
// sensitivity analysis restricted to states that passed a 
// law after 2008

* for different age groups and two different outcomes
foreach var of varlist pmu pru {
qui levelsof age, local(levels)
foreach a of local levels {
    di " "
    di "Future-passing states as control group"
    di "-> Age group = `: label (age) `a''"

	xtreg `var' i.year0 i.mmlf if age==`a' & treatex==1, i(state) fe cluster(state) cformat(%4.3f)
	outreg2 using mml-`var', sideway stats(coef ci) paren(ci) dec(2) cdec(1) noaster ///
		addtext(Year fixed effects, Yes, State fixed effects, Yes, Age group, `a')
		}
	}
		
log close analysis
exit
	
